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1.  Introduction 

This  work  is  concerned  with  the  development  of  an  efficient  parallel  Large 
Eddy  Simulation  (LES)  and  Detached  Eddy  Simulation  (DES)  capability 
using  unstructured  meshes.  The  advantages  of  unstructured  meshes  include 
flexible  modeling  of  complex  geometries,  adaptive  meshing  capabilities,  and 
homogeneous  data  structures  well  suited  for  massively  parallel  computer  ar- 
chitectures. On  the  other  hand,  unstructured  mesh  techniques  require  ad- 
ditional computer  resources  as  compared  to  cartesian  or  structured  mesh 
methods,  and  the  achievable  accuracy  of  the  particular  unstructured  mesh 
discretization  must  be  carefully  considered.  The  approach  developed  in  this 
work  is  based  on  an  existing  steady-state  unstructured  mesh  solver  which 
relies  on  agglomeration  multigrid  for  rapid  convergence  and  has  been  shown 
to  scale  well  on  inexpensive  personal  computer  (PC)  clusters  as  well  as 
on  massively  parallel  supercomputers  using  thousands  of  processors  [1].  A 
vertex-based  discretization  is  used,  where  the  flow  variables  are  stored  at 
the  vertices  of  the  mesh,  and  a single  edge-based  data  structure  capable 
of  handling  combinations  of  tetrahedra,  hexahedra,  prism  and  pyramids 
is  employed.  Spatial  discretization  is  achieved  through  a central-difference 
control-volume  formulation  with  scaled  matrix-based  artificial  dissipation 
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derived  from  an  upwind  Roe-Rieman  solver.  This  discretization  is  second- 
order  accurate  in  space.  The  baseline  steady-state  solver  is  extended  to 
an  unsteady  Reynolds-Averaged  Navier-Stokes  (URANS)  solver,  using  a 
second-order  accurate  three-point  backwards  difference  time  discretization 
[2].  At  each  physical  time-step,  the  parallel  agglomeration  multigrid  algo- 
rithm is  employed  to  drive  the  non-linear  (unsteady)  residual  to  conver- 
gence. The  agglomeration  algorithm  automatically  constructs  coarse  levels 
in  a pre-processing  phase,  based  on  the  graph  of  the  original  fine  grid. 
A preconditioned  multi-stage  iterative  scheme  is  then  used  on  each  grid 
level  to  drive  the  multigrid  algorithm.  Jacobi  (point-wise)  precondition- 
ing is  used  in  isotropic  regions  of  the  grid,  while  line  preconditioning  is 
used  in  highly  stretched  regions  of  the  grid  such  as  near  walls  where  high 
grid  stretching  is  required  to  resolve  thin  boundary  layers  [3].  The  steady 
and  unsteady  RANS  solver  employs  the  one  equation  turbulence  model  of 
Spalart  and  Allmaras  [4],  The  extension  from  RANS  to  LES  and  DES  is 
achieved  through  the  modification  of  the  length  scale  definition  d in  the 
Spalart- Allmaras  turbulence  model.  In  the  original  model,  d is  taken  as  the 
distance  at  a given  grid  point  to  the  closest  wall.  In  the  DES  extension 
proposed  by  Spalart  [5],  this  length  scale  is  replaced  by: 

dDES  = min(d,  Cdes , Ax)  (1) 

where  d represents  the  original  closest- wall  length  scale,  Cdes  represents 
a model  constant  taken  as  0.65,  and  Ax  represents  a measure  of  the  local 
grid  spacing.  For  unstructured  meshes,  Ax  is  taken  as  the  maximum  edge 
length  touching  a given  vertex.  In  regions  far  removed  from  walls,  this 
modification  to  the  turbulence  model  emulates  a simple  Smagorinsky  model 
for  LES,  while  reverting  in  a smooth  manner  to  the  well-established  Spalart- 
Allmaras  RANS  model  in  near  wall  boundary-layer  regions. 

The  resulting  solution  methodology  is  nominally  second-order  accurate 
in  space  and  in  time.  We  avoid  the  construction  of  higher-order  accurate 
spatial  operators  in  order  to  be  able  to  retain  all  the  discretization,  solu- 
tion, and  parallelization  techniques  previously  developed  in  the  steady-state 
solver  context.  However,  extra  care  must  be  taken  to  avoid  excessive  nu- 
merical diffusion  from  the  second-order  discretization  from  dominating  the 
turbulence  eddy  viscosity  levels  in  LES  regions.  With  this  in  mind,  the 
individual  convective  and  numerical  dissipative  terms  of  the  discretization 
are  evaluated  separately,  and  a (constant)  scaling  factor  is  applied  to  the 
dissipative  terms  which  provides  the  option  for  reducing  these  terms  from 
their  nominal  values  obtained  when  the  scheme  is  written  as  an  upwind 
Roe-Rieman  solver. 

The  agglomeration  multigrid  algorithm,  which  is  used  here  as  a non- 
linear implicit  solver  in  time,  removes  any  stability  restrictions  on  the  per- 
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missible  time-step  size,  thus  allowing  the  time-step  choice  to  be  determined 
solely  by  accuracy  considerations.  While  the  desirable  time-step  size  in  LES 
regions  may  be  small  enough  to  obviate  the  need  for  a fully  implicit  solver, 
in  the  DES  mode  very  high  Courant  numbers  will  almost  always  be  en- 
countered in  boundary  layer  regions,  where  the  solver  reverts  to  a RANS 
behavior,  thus  making  the  availability  of  an  implicit  solver  essential. 

In  the  following  sections,  the  solver  is  validated  first  in  unsteady  RANS 
mode  for  the  flow  over  a circular  cylinder.  Validation  in  the  LES  regime 
is  then  undertaken  by  computing  the  decay  of  isotropic  turbulence  in  a 
periodic  box.  Finally,  the  flow  over  a sphere  is  simulated  in  DES  mode  and 
compared  with  experiment  and  with  an  equivalent  simulation  using  the 
solver  in  an  unsteady  RANS  mode. 


2.  Unsteady  RANS  Flow  Over  a Circular  Cylinder 

The  flow  around  a circular  cylinder  is  a well-  known  case,  which  has  been 
widely  studied  computationally  and  experimentally.  This  case  is  used  as 
the  basis  for  validation  of  the  unsteady  RANS  solver,  and  for  assessing  grid 
resolution  and  time-step  requirements  for  accurately  predicting  the  vortex 
shedding  frequency  observed  in  the  cylinder  flow.  Two  different  meshes  of 
252,000  and  631,000  grid  points  and  three  different  time-steps  of  0.5,  0.25 
and  0.1  were  used.  The  time-step  is  non-dimensionalized  as  t — t0Uinf/d 
where  d is  diameter  of  the  cylinder  and  U{nf  is  the  freestream  velocity. 

The  one  equation  Spalart-Allinaras  turbulence  model  [4]  was  used  for  all 
calculations  in  fully  turbulent  mode.  In  all  cases  the  agglomeration  multi- 
grid strategy  was  used  with  four  levels.  The  Mach  number  is  0.2  and  the 
Reynolds  number  is  1200  for  this  case.  All  runs  were  performed  in  parallel 
using  16  processors  of  an  Intel  500  MHz  Pentium  III  PC  cluster. 

The  computational  domain  in  the  plane  normal  to  the  cylinder  span  has 
an  aspect  ratio  of  1 and  a side  length  of  100  cylinder  diameters.  A span  of 
two  cylinder  diameters  is  employed,  and  inviscid  (slip  velocity)  boundary 
conditions  are  applied  at  the  end-walls. 

Table  1 shows  the  Strouhal  Numbers  computed  for  each  mesh  and  each 
time-step  of  the  cylinder  flow  simulations.  Convergence  is  achieved  as  the 
time-step  is  reduced  and  the  mesh  size  increased.  A second-order  accurate 
convergence  behavior  is  observed  as  the  time-step  is  reduced,  validating  the 
accuracy  of  the  three-point  backwards  difference  scheme  used  to  discretize 
the  time-step.  From  the  smallest  time  step  results,  the  solution  can  be  seen 
to  be  grid  converged,  at  least  with  respect  to  the  prediction  of  the  vortex 
shedding  frequency.  The  computed  Strouhal  number  compares  very  well  to 
the  experimental  value  of  St  = 0.21. 
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3.  LES  Simulation  of  Isotropic  Decaying  Turbulence 

In  order  to  validate  the  solver  in  the  LES  mode,  the  simulation  of  decaying 
homogeneous  isotropic  turbulence  is  computed  and  the  energy  spectra  of 
the  flow  field  are  compared  with  experimental  results  from  Comte-Bellot 
and  Corrsin  [6].  Strelets  [7]  has  performed  similar  computations  which  were 
used  to  calibrate  the  value  of  the  Cdes  model  constant.  In  the  present  work, 
we  are  mainly  concerned  with  assessing  the  impact  of  numerical  dissipa- 
tion on  solution  accuracy,  since  the  scheme  is  only  second-order  accurate. 
Therefore,  the  value  of  Cdes  is  held  constant  at  0.65  for  all  computations, 
which  is  the  value  recommended  by  Strelets  [7].  Artificial  dissipation  levels 
are  varied  by  prescribing  different  values  of  the  dissipation  scaling  param- 
eter, and  by  varying  the  grid  resolution.  In  addition,  simulations  with  and 
without  the  turbulence  model  eddy  viscosity  are  performed  to  assess  the  ef- 
fect of  the  turbulence  model  on  overall  solution  accuracy,  and  to  determine 
whether  the  the  levels  of  eddy  viscosity  dominate  the  artificial  dissipation. 

The  computational  domain  consists  of  a cube  with  periodic  boundary 
conditions.  Because  no  wall  regions  are  present,  the  DES  model  operates  in 
the  LES  mode  throughout  the  domain  regardless  of  grid  resolution.  Compu- 
tations are  performed  on  a coarse  and  a fine  grid,  which  are  constructed  as 
unstructured  cartesian  grids  with  32  and  64  hexahedral  cells,  respectively, 
in  each  dimension  of  the  cubic  domain. 

In  order  to  compare  with  the  experimental  values  from  [6],  the  flow 
field  in  the  computational  domain  must  be  initialized  as  a solenoidal  field 
with  a prescribed  energy  spectrum  corresponding  to  that  measured  at  the 
initial  test  section  in  the  experiments  (up  to  the  cutoff  value  of  the  wave 
number  corresponding  to  the  grid  size).  The  initial  eddy  viscosity  field 
must  then  be  obtained  by  preconverging  the  turbulence  model  with  the 
flow-field  held  frozen.  Once  the  initialization  is  complete,  the  flow  field 
is  advanced  in  time  using  the  implicit  time-stepping  procedure  described 
above.  A time-step  of  0.01  is  employed,  where  time  is  non-dimensionalized 
as  t — t0(l.5u'2)1/'2 /L,  where  L is  the  box  length,  and  v!  represents  the 
initial  rms  average  velocity  fluctuation.  The  resulting  flow  fields  at  t = 0.87 
and  t=  2.0  are  postprocessed  to  obtain  the  energy  spectra,  which  are  then 
compared  with  corresponding  experimental  data. 

Figure  1 illustrates  the  computed  spectra  on  both  grids  at  two  time 
levels  for  the  nominal  value  of  the  artificial  dissipation  scaling  factor,  i.e. 
the  value  generally  employed  for  steady-state  calculations  in  RANS  mode. 
Clearly,  the  finer  scales  decay  more  rapidly  than  in  the  experimental  values. 
When  the  same  simulations  are  performed  with  the  eddy  viscosity  turned 
off,  little  difference  in  the  energy  spectra  is  observed,  suggesting  that  the 
eddy  viscosity  values  are  overwhelmed  by  the  levels  of  artificial  dissipation. 
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In  Figure  2,  the  same  computations  are  performed  with  a lower  scaling  of 
the  artificial  dissipation  terms  (0.25  of  the  nominal  value).  Substantially 
better  agreement  is  observed  at  all  scales,  although  the  finest  scales  are 
still  dissipated  faster  than  in  the  experimental  results.  However,  reasonably 
good  agreement  is  obtained  up  to  k=10  for  the  fine  grid,  and  both  grids 
agree  reasonably  well  at  lower  wave  numbers.  When  the  eddy  viscosity  is 
omitted  in  these  cases,  agreement  with  experiment  degrades,  particularly 
at  the  lower  wave  numbers  where  the  dissipation  is  lower  than  observed 
experimentally. 

4.  DES  Flow  Over  a Sphere 

The  current  solver  is  applied  in  DES  mode  to  predict  the  flow  around  a 
sphere.  A Mach  number  of  0.2  was  prescribed,  while  the  Reynolds  number 
was  set  to  104,  corresponding  to  the  sub-critical  regime  (laminar  boundary 
layer  separation),  which  is  similar  to  the  computations  performed  in  [8], 
which  provides  a comparative  basis  for  our  results.  An  unstructured  mesh 
of  767,000  vertices  is  employed,  in  a cubic  computational  domain  of  100 
sphere  diameters  in  total  length.  At  the  surface  of  the  sphere,  the  normal 
grid  spacing  is  10~4  sphere  diameters. 

Two  different  time-steps  of  0.1  and  0.05  are  used,  where  the  time-step 
is  non-dimensionalized  as  shown  in  equation  (1).  Four  multigrid  levels  are 
used,  and  each  physical  time-step  is  solved  with  25  multigrid  cycles,  result- 
ing in  a two  order  of  magnitude  reduction  in  the  residuals.  All  runs  were 
performed  on  a cluster  of  thirty-two  800  MHz  Pentium  III  PCs,  for  which 
the  convergence  of  a physical  time-step  could  be  achieved  in  approximately 
6 minutes  of  wall  clock  time. 

This  case  has  been  computed  in  both  DES  and  URANS  mode,  using  the 
Spalart-Allmaras  turbulence  model  without  modification  in  the  latter  case. 
The  time  history  of  the  drag  coefficient  is  shown  in  Figure  3 and  reveals  im- 
portant differences  between  URANS  and  DES.  The  mean  value  of  the  drag 
coefficient  in  both  cases  is  close  to  the  experimentally  reported  value  of  0.40 
from  Schlichting  [9].  However,  the  URANS  simulation  appears  to  damp  out 
most  of  the  oscillations  present  in  the  DES  run,  while  the  DES  runs  show 
a very  chaotic  oscillatory  pattern  quite  similar  to  the  solutions  obtained 
by  Constantinescu  et  al.  [8].  Spectral  analysis  of  the  time-dependent  drag 
coefficient  history  reveals  a peak  corresponding  to  a Strouhal  number  of 
0.1  which  is  not  in  agreement  with  the  values  0.18  to  0.2  reported  experi- 
mentally [8].  This  may  be  due  to  an  insufficiently  long  time  history  sample, 
since  less  than  three  full  periods  of  this  frequency  are  present  in  our  sample. 

In  Figure  4,  the  average  surface  pressure  computed  using  DES  is  seen 
to  provide  superior  agreement  with  experimental  results  at  Re=  165,000 
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reported  by  Achenbach  [10]  in  the  separated  region  over  the  URANS  re- 
sults, and  are  in  agreement  with  the  results  reported  by  Constantinescu  [8]. 
The  average  computed  separation  angle  of  81°  in  the  DES  case  compares 
reasonably  to  experimental  vale  of  82.5°  . 

5.  Conclusions  and  Further  Work 

This  work  prepresents  initial  efforts  at  developing  and  validating  a fully 
implicit  parallel  LES/DES  solver  based  on  unstructured  meshes.  In  the  near 
future,  we  intend  to  pursue  further  validation  studies  on  both  basic  flows 
using  finer  grids  and  time-steps,  and  more  complex  geometries  such  as  flow 
over  bluff  body  components  of  aerospace  vehicles.  Efficiency  improvement 
are  also  under  consideration  through  the  use  of  higher-order  time-stepping 
procedures  and  Krylov  acceleration  methods.  This  work  has  been  supported 
by  AFOSR  under  the  management  of  Dr.  Len  Sakell. 
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TABLE  1.  Computed  Strouhal  Number  for  Various  Grid  Sizes  and  Time  Steps 
for  RANS  Flow  over  Circular  Cylinder 


Grid  Size  (points) 

Time  Step  = 0.5 

Time  Step  = 0.25 

Time  Step  = 0.1 

0.252  million 

0.19249 

0.20304 

0.20833 

0.631  million 

0.19379 

0.20408 

0.20833 

vls2=20  ;Eddy  Vscoslty  ON 


vls2=20 ; Eddy  Viscosity  OFF 
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Figure  1.  Comparison  of  Computed  and  Measured  Energy  Spectra  for  Nominal  Artificial 
Viscosity  Levels  with  Eddy  Viscosity  (LEFT)  and  Without  Eddy  Viscosity  (RIGHT) 
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Figure  2.  Comparison  of  Computed  and  Measured  Energy  Spectra  for  Reduced  Artificial 
Viscosity  Levels  with  Eddy  Viscosity  (LEFT)  and  Without  Eddy  Viscosity  (RIGHT) 
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Figure  3.  Comparison  of  Computed  Drag  Coefficients  for  Sphere  using  DES  (LEFT) 
and  URANS  (RIGHT) 
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Mean  pressure  distribution  over  the  sphere. 
M-0.2,  Re  -10,000 


Figure  4 ■ Compaxison  of  Computed  Average  Surface  Pressure  Coefficient  using  DES 
and  URANS  versus  Experimental  Data 


